!--------------------------------------------!
!        A. Gazizov, LNGS, 27.12.2010        !
!  included in project ZN on 27.11.12, LNGS  !
!          coefficients on IR/UV of          !
!  $ \gamma + p \rightarrow e^+ e^- p $ and  !
!  $ \gamma + N \rightarrow N' + X $         !
!          divided by H_0                    !
!--------------------------------------------!

module pbdir
   use iruv, only: nEIR, EIR, WIR
!   use approx1x
   implicit none
   
   real(8),parameter::               &
        m_e = 5.1099892811d-4        &  ! mass of electron in GeV
      , m_p = 0.939d0                &  ! mass of proton in GeV
      , cH0 = 1.285d28               &  ! c/H_0 in cm
      , k_B = 8.617332478d-5         &  ! eV K^{-1}
      , T_0 = 2.72556d0	             &  ! T_{CMB} in K
	  , kT0 = K_B*T_0

   integer, parameter:: npair = 94   &
					  , nCMB  = 101	 &
                      , nNX   = 99   &
					  , nc    = 2    &  ! p->p or p->n; ...
					  , nm    = 6    &  ! [1,2..6]=(P,B,D,Bs,Ds,Dss)
                      , nz    = 51
   real(8), dimension(npair)    :: lgEeepir
   real(8), dimension(nNX  )    :: lgENXir, lnENXir
   real(8), dimension(nz   )    :: zzir
   real(8), dimension(npair,nz) :: beepir 
   real(8), dimension(2,npair,nz) :: bbseepir 
   real(8), dimension(nc,nm,nNX,nz):: PIR  ! pp/pn; 1,2..6; .. .. 

   real(8),dimension(nCMB):: lgEeepCMB, beepCMB  ! to check on CMB			   
   real(8),dimension(nEIR):: wtmp	! temporary array
						         
   real(8):: xiEg
   integer:: ic, im, ir_ready=0

!------------------
     contains
!------------------

!---------------------------------------------
Subroutine ReadBPDIR
   implicit none
   integer:: i, j
   open(1,file='Data/beepIR.dat')
      read(1,'(52E15.8)') (lgEeepir(i),      &
           (beepir(i,j),j=1,nz),i=1,npair)
   close(1)
   print*, " File beepIR.dat is read"
   do j = 1,nz
      zzir(j) = 1d-1*(j-1)
   enddo

   open(1,file='Data/PIRknz.dat')
      do i = 1,nNX
         read(1,'(E15.8)')	lgENXir(i)
         read(1,'(51E15.8)') (((PIR(ic,im,i,j),j=1,nz),im=1,nm),ic=1,nc)
      enddo
   close(1)
   print*, " File PIRkn.dat is read"

   open(1,file='Data/lnENX.dat')
      write(1,'(E15.8)') (lgENXir(i)*dlog(1d1),i=1,nNX)
   close(1)

   open(1,file='Data/PIRtkn.dat')
      do ic = 1,nc
         do im = 1,nm
		    do j = 1,nz
		       write(1,'(99E15.8)') (PIR(ic,im,i,j),i=1,nNX)
			enddo ! j
	     enddo ! im
	  enddo ! ic
   close(1)
   stop

!   open(1,file='Data/lnENX.dat')
!      read(1,'(99E15.8)')	(lgENXir(i),i=1,nNX)
!   close(1)



   ir_ready = 1

end Subroutine ReadBPDIR
!------------------------------

!------------------------------
Subroutine WrtCheck
   implicit none
   integer:: i
   open(1,file='Chck.dat')
      write(1,'(7E12.5)') (1d1**lgENXir(i)  &
	   , ((PIR(ic,im,i,1),im=1,3),ic=1,nc),i=1,nNX)

   close(1)
end Subroutine WrtCheck
!------------------------------

!---------------------------------------------------
! Get comoving \beta_{ee}(E,z)/H_0 on IR/UV Kneiske
! to be multiplied by (1+z)^3/Redshift(z1)
!
real(8) Function GetbeepIR(E,z)
   use msimsl, only: DQD2VL
   implicit none
   LOGICAL:: CHECK
   real(8),intent(in):: E, z
   real(8):: lgE, ztmp

   if(ir_ready.eq.0) call ReadBPDIR 
   if(z.gt.5d0) then
      ztmp=5d0
   else
      ztmp = z
   endif

   lgE = dlog10(E)
   if(lgE.lt.lgEeepir(1)) then
      GetbeepIR = 0d0
   else
      if(lgE.gt.lgEeepir(npair)) lgE = lgEeepir(npair)
      GetbeepIR = DQD2VL(lgE,ztmp,npair,       &
         lgEeepir,nz,zzir,beepir,npair,CHECK)
   endif

end Function GetbeepIR
!--------------------------------------------------

!--------------------------------------------------
Subroutine CalcBpair
   use approx1x, only: eep
!   use iruv, only: nEIR, EIR, WIR
   use gazlib, only: dgauss
   implicit none
             
   integer:: i, j, jj   
   real(8):: EE, R, xith  

   do i = 1,npair
      lgEeepir(i) = 1.37d1 + 1d-1*(i-1)
      EE   = 1d1**lgEeepir(i)           ! energy of proton in eV
      xiEg = 2d-18*EE/m_p               ! coef. to get Eg in GeV for a given $\xi$

      xith = m_p*eep(1)/(2d-18*EE)      ! eep(1)=1.0226d-3 is threshold in GeV of pair production  
!     xi_max = y_max = EIR(nEIR)/(1 eV)
      do j = 1, nz
         jj = (j-1)*4 + 1
         wtmp = WIR(:,jj)
         R = dgauss(Cpair,xith,EIR(nEIR),1d-8)
         beepir(i,j) = 2d-30*cH0*R
         print*, i, j, beepir(i,j)
      enddo   ! j
   enddo ! i
   open(1,file='beepIR.dat')
      write(1,'(52E15.8)') (lgEeepir(i),     &
           (beepir(i,j),j=1,nz),i=1,npair)
   close(1)
end Subroutine CalcBpair 
!--------------------------------------------------

!--------------------------------------------------
! xi = m_p*E_\gamma/(2*w_0*E)
real(8) Function Cpair(xi)
!   use iruv, only: nEIR, EIR
   use approx1x, only: eep, Avr1xpair
   use gazlib, only: divdif
   implicit none
   real(8),intent(in):: xi
   real(8):: Eg, avr, wint
   
   Eg = xiEg*xi    ! energy of nucleon in GeV
   if(Eg.le.eep(1)) then
      Cpair = 0d0
      return 
   endif

   call Avr1xpair(Eg,avr)
   if(xi.lt.EIR(1)) then
      wint = wtmp(1)
   else
      wint = divdif(wtmp,EIR,nEIR,xi,2)
   endif

   Cpair = xi*avr*wint

end Function Cpair
!---------------------------------------------------

!--------------------------------------------------
! Photonuclear interactions

Subroutine CalcPBD
   use approx1x, only: egg
   use gazlib, only: dgauss
!   use iruv, only: nEIR, EIR, WIR
   implicit none
             
   integer:: i, j, jj     ! 
   real(8):: EE, R, xith  ! 

   PIR = 0d0

! Energies of protons
   do i = 1,nNX  
      if(i.le.53) then 
	     lgENXir(i) = 15.8d0+0.05d0*(i-1)
      else
         lgENXir(i) = 18.5d0+0.1d0*(i-54)
	  endif
   enddo

   do i = 1,nNX
      EE   = 1d1**lgENXir(i)          ! energy of proton in eV
      xiEg = 2d-18*EE/m_p             ! coef. to get Eg in GeV for a given $\xi$
      xith = 1d18*m_p*egg(1)/(2d0*EE) !
	   
!     xi_max = y_max = EIR(nEIR)/(1 eV)
      do j = 1,nz
         jj = (j-1)*4 + 1
         wtmp = WIR(:,jj)
		 do ic=1,nc
		    do im = 1,3  ! calculate P, B, D
               R = dgauss(CPBD,xith,EIR(nEIR),1d-12)
               PIR(ic,im,i,j) = 2d-30*cH0*R
            enddo ! im
         enddo ! ic  
         write(*,'(2I3,6E12.5)') i, j, ((PIR(ic,im,i,j),im=1,3),ic=1,nc)
      enddo ! j
   enddo ! i
! Write result into file
   open(1,file='PIRkn.dat')
      do i = 1,nNX
         write(1,'(E15.8)')	lgENXir(i)
         write(1,'(51E15.8)') (((PIR(ic,im,i,j),j=1,nz),im=1,nm),ic=1,nc)
      enddo
   close(1)
end Subroutine CalcPBD 
!--------------------------------------------------

!--------------------------------------------------
! xi = m_p*E_\gamma/(2*w_0*E)
real(8) Function CPBD(xi)
   use approx1x, only: Avr1xMI
!   use iruv, only: nEIR, EIR
   use gazlib, only: divdif
   implicit none
   real(8),intent(in):: xi
   real(8):: Eg, avr, wint
   
   Eg = xiEg*xi                   ! energy of nucleon in GeV
   if(Eg.le.0.15d0) then
      CPBD = 0d0
	  return
   endif

   avr = Avr1xMI(Eg,im,ic)

   if(xi.lt.EIR(1)) then
      wint = wtmp(1)
   else
      wint = divdif(wtmp,EIR,nEIR,xi,2)
   endif

   CPBD = xi*avr*wint

end Function CPBD
!---------------------------------------------------

!--------------------------------------------------
! Recalculate \beta_{eep}(E) on CMB at z=0
Subroutine CalcBpairCMB
   use approx1x, only: eep
   use gazlib, only: dgauss
   implicit none
   integer:: i 
   real(8):: EE, R, xith  

   do i = 1,nCMB
      lgEeepCMB(i) = 1.6d1 + 5d-2*(i-1)	 ! test CMB
      EE   = 1d1**lgEeepCMB(i)           ! energy of proton in eV
      xiEg = 2d-18*kT0*EE/m_p            ! coef. to get Eg in GeV for a given $\xi$

      xith = m_p*eep(1)/(2d-18*kT0*EE)   ! eep(1)=1.0226d-3 GeV is the threshold of eep  
      if(xith.ge.3d2) then
	     R = 0d0
      else
         R = dgauss(CpairCMB,xith,3d2,1d-8)
      endif
      beepCMB(i) = 1d-30*cH0*(4.11d2/1.202055d0)*R
      print*, i, lgEeepCMB(i), beepCMB(i)
   enddo ! i
   open(1,file='beepCMB.dat')
      write(1,'(2E15.8)') (lgEeepCMB(i),beepCMB(i),i=1,nCMB)
   close(1)
end Subroutine CalcBpairCMB 
!--------------------------------------------------

!--------------------------------------------------
! xi = m_p*E_\gamma/(2*kT*E)
real(8) Function CpairCMB(xi)
   use approx1x, only: eep, Avr1xpair
   implicit none
   real(8),intent(in):: xi
   real(8):: Eg, avr
   
   Eg = xiEg*xi    ! energy of nucleon in GeV
   if(Eg.le.eep(1)) then
      CpairCMB = 0d0
      return 
   endif

   call Avr1xpair(Eg,avr)

   CpairCMB = -xi*avr*dlog(1d0-dexp(-xi))

end Function CpairCMB
!---------------------------------------------------

!---------------------------------------------------
Subroutine Derivs
   use gazlib, only: parab
   implicit none
   real(8),dimension(nNX):: tmp 
   real(8):: y,ys, yss
   integer:: i, j

   lnENXir = lgENXir*dlog(1d1)  ! arrary of ln(E)

   do ic = 1,nc
      do j = 1,nz
	     do i = 1, nNX
!            tmp(i) = dlog(PIR(ic,2,i,j))  ! for \beta'
            tmp(i) = dlog(PIR(ic,3,i,j))   ! for \delta' and \delta''
	     enddo ! i
		 do i = 1, nNX
            call parab(nNX,lnENXir,tmp,lnENXir(i),y,ys,yss)
!		    PIR(ic,4,i,j) = ys	! for \beta'
		    PIR(ic,5,i,j) = ys	! for \delta'
		    PIR(ic,6,i,j) = yss	! for \delta''
		    write(*,'(3(A,I2),E12.5)') ' ic=',ic, ' j=',j, ' i=',i, ys
		 enddo ! i 
      enddo	! j
   enddo ! ic

   open(11,file='Data/PIRknZ.dat')
!      write(11,'(99E15.8)')	(lgENXir(i),i=1,nNX)
!	  do j=1,nz
!		 write(11,'(99E15.8)')(((PIR(ic,im,i,j),i=1,nNX),im=1,nm),ic=1,nc)
!	  enddo ! j
      do i = 1,nNX
      write(11,'(E15.8)')	lgENXir(i)
      write(11,'(51E15.8)') (((PIR(ic,im,i,j),j=1,nz),im=1,nm),ic=1,nc)
      enddo	! i
   close(11)
   write(*,*) ' File Data/PIRknz.dat is rewritten'

end Subroutine Derivs 

!---------------------------------------------------


!---------------------------------------------------
Subroutine Deriveep
   use gazlib, only: parab
   implicit none
   real(8),dimension(npair):: lnEeep, tmp ! , BREAK, value
   real(8),dimension(npair,nz):: dbeepIR  ! derivatives d\ln\beta/d\lnE
   real(8):: y,ys, yss
   integer:: i, j

   lnEeep = lgEeepir*dlog(1d1)  ! arrary of ln(E)

   do j = 1,nz
      do i = 1, npair
         tmp(i) = dlog(beepir(i,j))
      enddo ! i
      do i = 1, npair
         call parab(npair,lnEeep,tmp,lnEeep(i),y,ys,yss)
		 dbeepIR(i,j) = ys
		 write(*,'(2(A,I2),E12.5)')  ' j=',j, ' i=',i, ys
      enddo ! i 
   enddo	! j

   open(11,file='Data/beepIR.dat')
      write(11,'(103E15.8)')	(lnEeeP(i),(beepIR(i,j),dbeepIR(i,j),j=1,nz),i=1,npair)
   close(11)
   write(*,*) ' File beepIRz.dat is written'

end Subroutine Deriveep

!---------------------------------------------------


end module pbdir